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I . INTRODUCTION 


Factorial  series  derived  from  the  Laplace  integral  converge  rapidly 
for  large  values  of  the  argument  and,  thus,  are  preferable  to  the 
corresponding  asymptotic  series.  However,  the  traditional  algorithm 
leads  to  very  large  numbers  and  must  be  modified  if  it  is  to  be  useful 
for  numerical  work.  One  procedure  for  scaling  the  large  Stirling 
numbers  which  occur  in  the  analysis  is  derived  below. 

Factorial  series  based  on  a Laplace  integral  evaluated  between 
finite  limits  will  generally  diverge,  so  that  an  alternate  procedure 
is  required.  One  method,  due  to  Hadamard,  is  to  expand  the  Laplace 
integral  in  a series  of  incomplete  gamma  functions.  The  resulting  series 
converge  rapidly  for  large  values  of  the  argument.  In  practice,  expan- 
sions in  terms  of  the  Kummer  function  are  more  convenient  for  computation. 
These  functions  are  closely  related  to  the  incomplete  gamma  function. 

Computer  programs  based  on  these  algorithms  will  be  used  to  check 
the  accuracy  of  the  BRL  subroutines  for  Bessel  functions  of  complex 
argument  and  integral  order.  This  is  necessary  as  tables  are  not  avail- 
able for  a sufficient  range  of  order  and  argument  to  make  a detailed  check 
by  comparison. 


II.  FACTORIAL  SERIES 

The  factorial  series  are  used  to  calculate  Kn(x) , Jn(x)  and  Y^Cx) 
K^fx)  can  be  expressed  in  terms  of  the  Whittaker  function  as1 


Kn(x) 


(fe) 


1/2 


W (2x)  , 
o,n  ’ 


where  the  asymptotic  expansion  for  the  Whittaker  function  is' 


W (2x)  = e 
o ,n 


-x 


i +£ 

m=l 


n2- C-l/2)2]  [n2-  (-3/2) 2] . 

. . rn2- (1/2  - m)2]  ) 

m! (2x)m 

1 

1 Handbook  of  Mathematical  Funotionsy  NBS55,  U.S.  Government  Printing 

Office,  1964,  p.  377. 


2 Modern  Analysis,  E. 
Cambridge,  England, 


J.  Whittaker  and  G. 
1927,  p.  343. 


N. 


Watson, 


University  Press, 
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This  asymptotic  expansion  was  derived  from  a Laplace  integral  evaluated 
between  zero  and  infinity  and  involves  only  negative  integral  powers  of 
the  argument. 


For  n = 0, 


K0(x)  = 


' -X  j l2  12~32  12,32,52 

1!(8x)  2 ! (8x) 2 3!(8x)3 


For  n = 1 , 


1/2  k A. 

' e-  E 4 

j =o  xJ 


2 2 2 

-x  . 1 * 3 r -3’5  1 ‘3  5*7 

1!(8x)  2 ! (8x) 2 3!(8x)3 


r i _M1/2  -x  j,  1 
K1W  \2x/  e j1  + l!(i 

/ \l/2  k B. 

= {k)  £ -i 

V f ]=o  xJ 


A computer  tabulation  of  the  first  fifty  of  these  coefficients  is 
shown  in  Table  I. 

These  series  can  be  summed  by  convergent  factorial  series  using  an 
algorithm  described  by  Wasow: 


-p  _ v r-p+1 


- E 

r=p- 1 


x(x+l)(x+2)  . . . (x+r) 


where  r denotes  the  Stirling  numbers  of  the  first  kind.' 


'•  K<,w  "(h)1 


e~X  S 


where  S = 1 - Tr7o~T  + T 
o 1 ! (8x)  o 


k A.  A„  A 


IN.  r\  . Ha  rv  — 

= T 4 + 4+  . . . 


j=2  xJ  x x 


Asymptotic  Expansions  for  Ordinary  Differential  Equations , W.  Wasow , 
Intersoienoe  Publishers , John  Wiley , NY,  1965,  p.  330. 
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0 

0. 100000000000000E 

01 

1 

-0. 125000000000000E 

00 

2 

0.703125C00000000E- 

-01 

3 

-0.732421875C00000E- 

-01 

4 

0. 1121520996G9375E 

00 

5 

-0.227108C01708984E 

00 

6 

0.5725C1420974731E 

00 

7 

-0. 172772750258446E 

01 

8 

0.607404200127348E 

01 

9 

-0. 243805296995561 E 

02 

10 

0. 110017140269247E 

03 

11 

-0. 551335896122021E 

03 

12 

0.303 80905 1092 238 E 

04 

13 

-0. 182577554742932E 

05 

14 

0. 118838426256783E 

06 

15 

-0.832859304016289E 

06 

16 

0.625 2951 49343480 E 

07 

17 

-0. 500695895319889E 

08 

18 

0.425939216504767E 

09 

19 

-0. 383625518023043E 

10 

2C 

0. 364684008070656E 

11 

21 

-0. 364901081884983E 

12 

22 

0. 3833534661 39 394E 

13 

23 

-0.421897157028410E 

14 

24 

0.48540 146 86 85 290 E 

15 

25 

-0.582724463156691E 

16 

26 

0.728685734937766E 

17 

27 

-0.947628809926011E 

18 

28 

0. 127972194197597E 

20 

29 

-0. 179216232305170E 

21 

30 

0.259938210272623E 

22 

31 

-0. 390012129203400E 

23 

32 

0.604671 14 8 7 53 240 E 

24 

33 

— 0.96 77028801 06 985E 

25 

34 

0. 1597C6552529421E 

27 

35 

-0.271558177354491E 

2« 

36 

0 . 475  321 10 1404 162E 

29 

37 

-0. 8557385639R0669E 

30 

38 

0.158339783631292E 

32 

39 

-0.300896338830105E 

33 

4C 

0.5P6841890824589E 

34 

41 

-0.  U7386269685980E 

36 

42 

0.240676789246046E 

37 

43 

-0.505491221599616E 

38 

44 

0. 108694973189986E 

40 

45 

-0.239159134066077E 

41 

46 

0. 538173040543799E 

42 

47 

-0. 123794112437854E 

44 

48 

0.290948402279072E 

45 

49 

- 0.698 3 503 870 00965 E 

46 

0.100000000000000E  01 
0.375000000000000E  00 
-0.117187500000000E  00 
0.10253 9062500000E  00 
-0.144195556640625E  00 
0.277576446533203E  00 
-0.676592588 42 4683E  00 
0.199353173375130E  01 
-0. 6883914268 10995E  01 
0.27248827311 268 5E  02 
-0.121597891876536E  03 
0. 6038440 7 6705070E  03 
-0.3 30227229 44 8085E  04 
0.197183759122366E  05 
— 0.127641272646175E  06 
0.890297 87 67070 6 8E  06 
-0. 665636771881769E  07 
0.531041101096852E  08 
-0.450278600 30 5039E  09 
0.404362032 510 775E  10 
-0.3833857520742 7 9E  11 
0,3827011 3465986 IE  12 
— 0.401183859913320E  13 
0. 44064 P141785228E  14 
-0. 50605685033 1473E  15 
0, 606509 135 122 2 70 E 16 
-0.757261646111796E  17 
0.983388387659068E  18 
-0. 132625728532056E  20 
0.185504521157983E  21 
-0.26874 96 75027628E  22 
0.40279 94 12128102E  23 
-0.62386 70 58237470E  24 
0.997478353341046E  25 
-0. 164473912306419E  27 
0.27942 94 288720 12E  28 
-0.488710428204280E  29 
0. 8 791 83456 1445 2 3E  30 
-0.162562177861459E  32 
0.308711828150368E  33 
-0.60169«647554326E  34 
0.120284696 09 7979E  36 
-0. 24647 6229950770E  37 
0 . 5 17385 1 32696078 E 38 
-0.111193708205847E  40 
0. 2445 3 3 4966 293 59 E 41 
-0. 550001019456850E  42 
0.126456351415012E  44 
-0.297073631R00736E  45 
0 . 71 274936405 2532E  46 


Table  I.  Coefficients  for  Asymptotic  Series 
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3 A3 


(x(x+ 


1)  + x(x+l) (x+2)  x(x+l) (x+2) (x+3) 


Therefore,  Tq  can  be  expressed  as 


t = y 

o ^ x(x+l)  . . , (x+r)  * 


where  V = A_  T . + A_  T _ + A.  T _+ 
o,r  2 r-1  3 r-2  4 r-3 


These  coefficients  can  be  calculated  and  stored  in  the  memory  of  the 
computer  for  recall  on  demand.  The  calculations  for  these  coefficients, 
involving  Stirling  numbers,  lead  to  very  large  numbers  in  the  computa- 
tion of  high-order  terms. 


Since  the  Stirling  numbers  are  always  greater  than  or  equal  to  one, 
we  modified  them  for  optimal  use  of  the  full  range  of  the  computer. 


The  Stirling  numbers  were  modified  in  the  following  way: 


S^  = F I^/Cv-l) ! , F = scale  factor,  such  as  10*'5 

S1  = F 
o 

so  = 

sl  = sl\  ♦sT1/<v-1> 

The  scale  factor  and  the  number  of  modified  Stirling  numbers  which  can 
be  calculated  are  machine-dependent.  The  computers  at  BRL  have  a range 

from  lO*^5  to  10155,  single  precision,  which  is  larger  than  the  range 

125 

of  most  computers.  As  can  be  seen  from  Table  II,  for  F = 10  and 

n = 150,  the  modified  Stirling  numbers  range  from  10  135  to  lO1*-5.  The 
process  of  scaling  the  Stirling  numbers  in  this  way  must  then  be  reversed 
in  calculating  each  term  of  the  factorial  series. 

By  this  transformation  we  obtained  accurate  results  (15  significant 
digits)  for  x ^ 6 by  summing  150  terms.  Similar  accuracy  could  be  ob- 
tained on  most  computers  using  double  precision. 
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1 0.262541431038901-1 

4 0.594416307877133-1 

7 0.620095188981095-1 

10  0.141690995161720-1 

13  0. 113521536148685E- 

16  0.4031C5100240183E- 

19  0.730187886089758E- 

22  0.740917C60108629E- 

25  0.450102130711339E- 

28  0. 172006626264944E- 

31  0.429520654391673E- 

34  0. 72214995206318 0E- 

37  0.837326275585929E- 

40  0.6827278C6061237E- 

43  0.397758943581573E- 

46  0. 167774C15729906E- 

49  0. 5179372 10655300 E- 

52  0. 118075771165955E 

55  0.200240793190783E 

58  0.2 54 1074421 06 383E 

61  0.242426888397106E 

64  0. 174495356761975E 

67  C.950002380342511E 

70  0.391821615456175E 

73  0.1225C8995716187E 

76  0.290319945584933E 

79  0.620899483162194E 

82  0. 7062546321 17985E 

85  0 .721564602846948E 

88  0. 55343P118324812E 

91  0. 3171672946794235 

94  0. 135025685209439E 

97  0.424C6C880287624E 

ICO  0.974387656698C37E 

1C3  0.162215177246175E 

1C6  0. 193432814981343E 

1 C 9 0. 1630C54839C6049E 

112  0.955495084811848+1 

115  0.3823646C8112042+1 

118  0.102158281160150+1 

121  0.197415994653956+1 

124  0.193865171243942+1 

127  0.128034645087299+1 

130  0.485782753348748+1 

133  0.991494859813093+1 

136  0.996598263274595+1 

139  0.435467363384560+1 

142  0.6844C32487P7668+1 

145  0.279 759712120229+1 

148  0.147742753080902+1 


35  0.293390049185972-131 

24  0.161631807256184-120 

14  0.937263899085794-111 

04  0.145745824899596-101 

95  0.872724632875019E-93 

87  0.244485921578491 E-84 

79  0.361878183982837E-76 

71  0.3075 06122625604E-68 

63  0. 1592 90231 850353 E-60 

55  0.526244225186191E-53 

48  0. 114831115329136E-45 

41  0.1701611725 44539E-38 

34  0.175105193623330E-31 

27  0.127434483546970E-24 

20  0. 66576645797338  IE-18 

13  0.252791612712235E-11 

07  0.704743584443506E-05 

00  0.145466542291410E  02 

06  0.223832004301461 E 08 

12  0.258158475 955319E  14 

18  0.2241 374358172 19E  20 

24  0.146958814419687E  26 

29  0.72926885372 65 16E  31 

35  0.274247724384982E  37 

41  0.781 85678301 7240E  42 

46  0.1 68899572 458664E  48 

51  0.2760965230914755  53 

56  0. 3407683 70204050 E 58 

61  0.31656R545939851E  63 

66  0.220455227036434E  68 

71  0.1 14501970441310E  73 

76  0.440823733512481E  77 

80  0.1 24873374341 855E  82 

84  0.258006379428166E  86 

89  0.3 84 83 63 75 802 625 E 90 

93  0.409407659018037E  94 

97  0.306267211155246E  98 

00  0.158431163016062+102 

04  0.555606274242301+105 

08  0.129004865734295+109 

11  0.192739054938949+112 

14  0.178944908244716+115 

17  0.988621557189148+117 

19  0.307581787025027+120 

21  0.501359095737511+122 

23  0.388006214031004+124 

25  0.123698885459199+126 

26  0.129991588080456+127 

27  0.300552651141317+127 

27  0.558451392197723+126 


0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 

0. 


Table  II.  Modified  Stirling  Numbers  for  n = ISO 
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162469629570885-127 
348403288776085-117 
122805989253782-107 
13499468151 1499E-98 
61 69648 62448953 E- 90 
138176751796 042 E- 81 
16865 128328499  IE-7  3 
12081072501 1220E-65 
536288624778063 E- 58 
1537591 7819021 5E-50 
29408952486491 6 E-4  3 
38 5045773882 69 4E-36 
352369350644976E-29 
2292666345940 32 E-2 2 
10755538788 3677E-1 5 
3680433898311 36E-09 
927 44 16 56 806729 E-03 
17 3458 45 3446 4 50 E 04 
242317917476252E  10 
25412 95 80523208E  16 
2008639 1046642 IE  22 
119996196936958E  28 
542840153293270E  33 
1861389 69 273117E  39 
483840847139701E  44 
95264530631 1109E  49 
141844276858788E  55 
159312784381071E  60 
134510966747108E  65 
850010949437138E  69 
399846 12 1 274463E  74 
13909552299756 IE  79 
355050280369826E  83 
6588880381 38 901E  87 
8793483 40 6 49 532E  91 
8 332934 16906864E  95 
552343491325650E  99 
251599210263367+103 
771214486559342+106 
155126469001113+110 
198618366362429+113 
155930607458550+116 
716280594664574+118 
181290891228498+121 
233459123119021+123 
136972391170442+125 
311019186727229+126 
209410799774947+127 
248534593164330+127 
100000000000000+126 


Similarly,  K^x)  = e A si 


= 1 + 

1*3 

1!  (8x) 

oo 

V 

- La 

r=l 

x(x+l) 

+ T, 


x+1)  . . . (x+r)  * 


where  = B^^  * B^^  * * . . . 

The  results  for  K^x)  were  equally  accurate. 

.The  asymptotic  series  for  the  ordinary  Bessel  functions,  x <_  25, 


are : 


II 

'rt 

*> J 

o 

(2  V/2 

W 

[Pq(x)  cos(x  - 

•JjM  = I 

f 2 v/2 

X/ 

[P  (x)  cos(x  - 

Y0(x)  = | 

(fT 

[Pq(x)  sin(x  - 

V*)  = 1 

[Pj(x)  sin(x  - 

rv^n.  1 _ . 

2 2 
1*3 

12’32*52*72 
+ . 

r)l 


2 ! (8x) 2 4 ! (8x)4 


k C. 

= £ 


j=o  X2j 

, , l2  12,32'52 


l2-32,52’72‘92 
5!  (8x)5 


k 

- £ 

j=o  X 


D. 

_J 

2j  + l 


4 Bessel  Functions,  Part  I , published  by  British  Association  for  the 
Advancement  of  Science , University  Press}  Cambridge } England. , IPS?. 

p.  202. 
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F 


Note  that  Cq  = | Aq | , Cx  = -|A2|,  . . . , C.  = (-1)3 |A  | 
and  Dq  = -lAjI,  Dx  = |A3|,  . . . , D.  - (-1) 3 + 1 |A2 j + J | 


Similarly, 


J=0  X J 


F. 


and  a.  E -Jrjy 


J=0  X 


And,  again,  Eq  = I B0 I , Ej  = - | B2 | , . . . , = (-l)3|B2j| 


F0  ■ |Bj|.  Fj  - , F.  - (-DhB^I 


For  the  ordinary  Bessel  functions,  x > 25, 
J (x)  = G(x)  sin(x)  + H(x)  cos(x) 

(x)  = M(x)  sin(x)  - N(x)  cos(x) 

Yo(x)  = H(x)  sin(x)  - G(x)  cos(x) 

Y (x)  = -N(x)  sin(x)  - M(x)  cos(x)  , 
where  G(x)  = (ttx)~1</2  [PQ(x) -Q0(x)  ] 

H(x)  = (itx)'1^2  [Pq(x)+Qo(x)  ] 

M(x)  = (ttx)_1//“  [P-,  (x)+Qj  Cx)] 

N (x)  = (ttx)'1/2  [Pj  (x) -Qj  (x) ] 


So,  for  x > 25,  the  same  coefficients  are  merely  arranged  in  a different 
manner. 

As  before,  the  results  obtained  were  accurate  to  15  significant 
digits  for  x > 6 by  summing  150  terms.  A sample  tabulation  of  the 
ordinary  Bessel  functions  from  the  computer  is  shown  in  Table  III. 

We  attempted  to  calculate  I (x)  in  the  same  manner,  but  the  factorial 
series  diverged. 

III.  HADAMARD  SERIES 

The  factorial  series  for  calculating  IQ(x)  and  Ij(x)  diverge 
since  the  Laplace  integrals  representing  these  functions  are  taken 
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1422*47282678085-01  0. 20510403861 3523E  00  0.2054642960389186  00  0.2107362803687356-01  FACT. SER I ES 

1422447282678085-01  0.205104038613522E  00  0 .205464296038919E  00  0.210736280368736E-01  SUBROUTINE 


j 


, 


: 

\ 


5 


t i 


r • 

f>  1 

f : 

f 

1 

L 


between  finite  limits  and,  therefore,  cannot  be  expanded  according  to 
the  previous  algorithm.  The  Hadamard  series,  useful  for  large  x,  was 
used  instead  and  has  been  programmed. 

In(x)  can  be  expressed  by:^ 


j (x)  = (x/?) 

nW  T(n+l/2)r(l/2) 


P 

/ x cos0  . 2n 
/ e sm 

Jo 


After  expansion  and  term-by-term  integration,  the  Hadamard  series  can 


(l/2-n)m  y(n+m+l/2,  2x) 
m! (2x)m 

where  y denotes  the  incomplete  gamma  function  and  (1/2-n)  denotes 
Pochhammer's  symbol. 

(a)n  = a(a+l) (a+2)  . . . (a+n-1),  (a)Q  = 1 

Each  term  in  the  expansion  of  these  series  contains  the  incomplete  gamma 
function,  which  is  expressed  below  in  terms  of  the  Kummer  function. ° 

y(a,x)  = a * xa  e X M(l,  1+a,  x)  , 
where  M denotes  the  Kummer  function. 

Hence,  after  substituting  and  simplifying,  we  have 

P'xf?Y'in  ^ C1/2- h)—  M (1,  n+m+3/2,  2x) 

I fxl  = - J > 

nl } r(n+l/2)T(l/2)  (n+m+1/2)  m! 


then  be  written  in  the  form 

, (x)  . (2x)-^2  t 

nl  J r(n+l/2)r(l/2)  m=o 


The  solution  of  these  series  is  straightforward  and  presented  no 
problems  in  overflowing  the  memory  of  the  computer.  The  calculation 
of  the  Kummer  function  required  many  terms  (250  terms  for  x=75)  to  get 
the  required  accuracy.  The  solutions  of  the  Hadamard  series  seem  to 
have  the  correct  convergent  behavior.  A sample  computer  tabulation  is 
shown  in  Table  IV. 


Theory  of  Bessel  Functions , 2nd  Ed.  > G.  N.  Watson,  Macmillan  Co., 
N.Y.,  1948,  p.  204. 

Handbook  of  Mathematical  Functions , NBS  55,  II.  S.  Government  Printing 
Office,  1964,  pp.  262,  504. 
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Table  IV.  Computer  Tabulation  of  I (x)  and 


The  results  were  good  but  not  as  accurate  for  moderate  argument  as 
we  had  hoped.  We  obtained  15  significant  digits  for  x ^ 17  by  summing 
25  terms  or  less  in  the  Hadamard  series.  This  is  not  much  better  than 
the  asymptotic  series  given  by  * 


>0w 

Ij(x) 


x .2  2.  2 2.  2.-2 

1 + T^  + W + 7W 

x 2 2 2 

e 1 ' 3 lz-3’5  r-3  '5*7 

- 1/2  " 1 ! (8x)  ' .2  ' .3 

(2itx)  j 2 ! (8x)  3!  (8x) 


When  the  asymptotic  series  were  programmed,  we  obtained  15  significant 
digits  for  x _>  19.  However,  the  Hadamard  series  does  provide  an 
independent  check  on  the  accuracy  of  the  asymptotic  series  used  in 
our  Bessel  function  subroutine. 


P 


ir;  « 

I i 


IV.  DISCUSSION  AND  CONCLUSIONS 

Factorial  series  are  an  effective  method  for  calculating  modified 
Bessel  functions  of  the  second  kind  and  related  functions.  Calculations 
have  been  limited  to  real  arguments  in  this  report;  however,  it  is  antici- 
pated that  extension  of  the  algorithm  to  complex  argument  will  not 
present  any  major  difficulties. 

On  the  other  hand,  the  Hadamard  series  does  not  present  any  advan- 
tages over  the  usual  asymptotic  series,  and,  consequently,  extending 
the  algorithm  to  complex  arguments  is  not  planned.  An  expansion  of 
the  Laplace  integrals  for  I (x)  and  I ^ (x)  in  terms  of  the  incomplete 
beta  function  is  now  being  developed  and  should  overcome  difficulties 
encountered  with  the  Hadamard  series. 
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